Extinction of an infectious disease: a large fluctuation in a non-equilibrium system 
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We develop a theory of first passage processes in stochastic non- equilibrium systems of birth- 
death type using two closely related epidemiological models as examples. Our method employs the 
probability generating function technique in conjunction with the eikonal approximation. In this 
way the problem is reduced to finding the optimal path to extinction: a heteroclinic trajectory of 
an efi'ective multi-dimensional classical Hamiltonian system. We compute this trajectory and mean 
extinction time of the disease numerically and uncover a non-monotone, spiral path to extinction 
of a disease. We also obtain analytical results close to a bifurcation point, where the problem is 
described by a Hamiltonian previously identified in one-species population models. 
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Statistics of large fluctuations in stochastic non- 
equilibrium systems has received much recent attention 
While the equilibrium fluctuation probability is de- 
termined by the Boltzmann distribution, there is no sim- 
ilar general principle away from equilibrium. The un- 
derlying reason is the absence of time reversal symme- 
try between the relaxation and excitation dynamics in 
out-of-equilibrium systems. As a consequence, the most 
probable fluctuation path is not determined by the relax- 
ation trajectory of the underlying deterministic system. 

An important class of stochastic non-equilibrium sys- 
tems is reaction kinetics, or birth-death, systems 
Rather than being caused by external factors, the noise 
in these systems is intrinsic, as it originates from dis- 
creteness of the reacting agents and random character of 
their interactions. When the typical number of agents 
is large, the Fokker-Planck (FP) approximation to the 
master equation, see e.g. Ref. [3], can accurately de- 
scribe small deviations from the probability distribution 
maxima. It fails, however, in determining the probability 
of large fluctuations 0, Therefore, developing ade- 

quate theoretical tools for dealing with large fluctuations 
is an important task. 

One of the areas where the birth-death models have 
been very successful is mathematical epidemiology, see 
Refs. [6[. In this Letter we consider two closely re- 
lated models of spread of disease in a population. Al- 
though they have served as standard multi-population 
epidemiological models, the analysis of large fluctuations 
in each of these models has not been satisfactory. We will 
use the two models as prototypical examples of multi- 
dimensional stochastic non-equilibrium systems. 

Observing the dynamics of a disease in a finite popula- 
tion, one notices the remarkable phenomenon of extinc- 
tion of the disease in a finite time. The expected time 
to extinction (and the possibility to affect it) is of great 
practical interest. Here we develop an efficient theoretical 
approach capable of computing, among other things, this 
quantity. The approach employs the probability gener- 
ating function formalism in conjunction with the eikonal 



approximation. In this way the problem is reduced to the 
dynamics of an effective classical Hamiltonian system. 
The intrinsic-noise-induced extinction of the disease pro- 
ceeds, with a high probability, along the optimal path: 
a special (heteroclinic) trajectory in the phase space of 
the classical Hamiltonian flow. An additional challenge 
of this type of problems is in the fact that the emerging 
multi-dimensional Hamiltonian flows are generally non- 
integrable. We compute the optimal path, and the mean 
extinction time of the disease, numerically and also ob- 
tain analytical results close to a bifurcation point. 

Model. Let us consider two stochastic epidemiologi- 
cal models: the endemic SI model and the endemic SIR 
model. In the SI model the host population is divided 
into two dynamic sub-populations: Susceptible (S) and 
Infected (I). The model is specified by the set of reac- 
tions and their rates given in Table 1. We can always 
represent the renewal rate (an independent parameter of 
the model) as fj,N, where N scales as the total popula- 
tion size in a steady state. Taking ^/ > /i, one allows an 
increased death rate of the infected. 



Event 


Type of transition 


Rate 


Infection 
Renewal of susceptible 
Death of susceptible 
Death of infected 


S^S-1, I^I + l 
S + 1 
S-1 
I ^ I -1 


{P/N)SI 
liN 
fiS 
fill 



TABLE I: Transition rates for the stochastic SI model 



The endemic SIR model deals, in addition to the S- 
and /-sub-populations, with a third sub-population: Re- 
covered (R) , with the recovery rate 7/. It is assumed that 
the recovered cannot become susceptible. The death rate 
of the recovered is iirR. The endemic SIR model (which 
generalizes the original SIR model: the one without re- 
newal and death) gives a satisfactory description to the 
spread of measles, mumps and rubella [^. 

Let us briefly review the deterministic, or mean-field 



2 



version of the SIR model: 



S ^ fiN - fiS - {f3/N)S I , 
i = -^hI -iI+{(3/N)SI , 
R = -^irR + jI. 



(1) 
(2) 
(3) 



As the dynamics of S and / decouples from that of R, 
the SIR model is effectively two-population, and we will 
not deal with the i?-dynamics. Furthermore, one imme- 
diately notices that, by putting /i/ -t- 7 = F, the S- and 
I-dynamics in the SIR model becomes identical to that 
in the SI model, up to interchange of /x/ and F. This also 
holds for the stochastic versions of the two models, and 
so we can treat them on equal footing, using F for the 
effective death rate constant of infected. 

For a sufficiently high infection rate, /3 > F, there is 
an attracting fixed point 



A^(/3-r) 
I3r 



N 



(4) 



which describes an endemic infection level, and an un- 
stable fixed point S — N, I — which describes an unin- 
fected population. At ^ < 4 (/3-F)(F//3)^ the attracting 
fixed point is a stable focus, while in the opposite case 
it is a stable node. The inverse of the real part of the 
eigenvalues (for the focus), or the inverse of the smaller 
of the eigenvalues (for the node) yields the characteristic 
relaxation time Tr towards the "endemic point" . 

The stochastic formulation of the SI and SIR mod- 
els accounts for the demographic stochasticity and ran- 
dom character of contacts between the susceptible and in- 
fected. The master equation for the probability Pn,m(t) 
of finding n susceptible and m infected individuals has 
the form 

Rn,m [A^(-^n — l,m Rji.m) ^" (^ l)-f7i+l,m '^Rn,m\ 

+ F [(771 + 1)P„ 

,m+l '^-Pn,m\ 

+ {l3/N)[{n + l){m-l)Pn+i ] ,(5) 

and the total population size is fluctuating in time. We 
will be interested in the regime where the fluctuations are 
relatively weak. In this case, after the relaxation time 
Tr a long-lived (quasi-stationary) distribution is formed 
that has a bi-variate gaussian peak with relative width 
^ around the stable state ^ of the mean-field de- 

scription 0, H, Q . The long-time behavior of the stochas- 
tic model is quite remarkable: due to a rare sequence of 
discrete events the disease goes extinct in a finite time. 
Given that a major outbreak of the disease occurred, 
what is the mean extinction time r of the disease? For 
the endemic SIR model this question was addressed pre- 
viously 0, @] the framework of the van Kampen sys- 
tem size expansion that brings about the approximate 
FP equation 0. Our approach considerably (exponen- 
tially) improves on these earlier results. In the regime we 
are interested in r is exponentially large compared with 



the relaxation time Tr- The presence of the large pa- 
rameter facilitates the use of the eikonal approximation: 
either directly in the master equation, as suggested by 
Dykman et al. To*], or in the evolution equation for the 
probability generating function, as suggested by Elgart 
and Kamenev 

Probability generation function and eikonal approx- 
imation. We adopt the latter approach and intro- 
duce the probability generating function G(ps,pi,t) = 
YlZm=oP'sPTPn.rn{t)- Ouce G{ps,Pi,t) is fouud, the 
probabilities Pn,m{t) are given by the coefficients of its 
Taylor expansion around ps = Pi = 0- Using the mas- 
ter equation ([5|), we obtain an evolution equation for G: 
dtG = HG with the effective Hamiltonian operator 



/i(ps-l)(iV-5p,)-F(p,-l)9j 



{PlN){ps-pi)pidl 



(6) 



In contrast to the FP equation this equation is exact [ll| . 

The eikonal ansatz is G{ps,pi,t) ~ cxp[—S{ps,pi,t)], 
where 5^1. Neglecting the second derivatives of S 
with respect to ps and pj, we arrive at a Hamilton- Jacobi 
equation dtS + H = in the p-space with the classical 
Hamiltonian H{S, I,ps,pi)' 

H = ^,(ps-l){N-S)-Tipi-l)I-iP/N)ips 



-pi)piSI , 
(7) 

where S — —dp^S and / — —dp^S. The structure of four- 
dimensional (4D) phase space, defined by the Hamilto- 
nian (O, provides a fascinating and instructive insight 
into the disease extinction dynamics. As H does not 
depend explicitly on time, it is an integral of motion: 
H{S, I,ps,pi) — E — const. All the mean-field trajecto- 
ries, described by Eqs. (HJ and lie in the zero-energy, 
E = Q, two-dimensional plane ps = pi — 1. The at- 
tracting fixed point (|4|) of the mean-field theory becomes 
a hyperbolic point A = [S,I,1, 1] in the 4D phase space 
with two stable and two unstable eigenvalues (the sum of 
which is zero) and respective eigenvectors. There are two 
more zero-energy fixed points in the system: the point 
C — [N, 0, 1, 1] which is present in the mean-field de- 
scription and the non-mean-field point B = [N, 0, 1, F//?] 
which we call fiuctuational. Both of them are hyper- 
bolic and describe extinction of the disease. The pres- 
ence of a fiuctuational fixed point, related to extinction, 
is characteristic of a class of stochastic birth-death sys- 
tems 0,0, Em. 

The most probable sequence of discrete events, bring- 
ing the system from the endemic state to extinction of the 
disease, is given by the optimal path that minimizes the 



action 5 10|)llJ|- The optimal path must be a zero-energy 



heteroclinic trajectory. This trajectory exits, at t = —00, 
the "endemic" point A along its two-dimensional unsta- 
ble manifold and enters, at t = cxd, the fiuctuational dis- 
ease extinction point B, along its two-dimensional sta- 
ble manifold. As in one-dimensional birth-death systems 
O'^^ '^^^ show that there is no trajectory going 
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directly from A to C. Therefore, the fluctuational ex- 
tinction point B, not present in the mean- field dynamics, 
plays a crucial role in the disease extinction. 

Up to a pre-exponent, the mean extinction time of the 
disease is t oc Trexp(<So) where 



{psS + pii)dt, 



(8) 



and the integration is performed along the (zero-energy) 
optimal path. As in any generic multi-dimensional 
Hamiltonian system, the optimal path can be computed 
only numerically. In the following we present two typi- 
cal examples of such computation, and also consider an 
important limit when the computation can be performed 
analytically, by exploiting time scale separation. First we 
introduce new coordinates x = S/N — 1 and y — I/N, 
time t = fit^ momenta px.y — psj — 1 and bifurcation 
parameter S— 1~T/(3,0<S< 1. The action ([8|) can 
now be rewritten as Sq — Na, where (t{K, 6) is the action 
along the optimal path, generated by the Hamiltonian 

H = -P,X ~K[{1~ S)Py + [p., ~ Py){Py + l) {x + l)] ^ 

(9) 

and K = j3/ ji > 1. The fixed points A, B, and C become 
^-^,0,0 , [0,0,0,-5], and [0,0,0,0], 



-s, 



K{1 
respectively. 

Optimal path and action: numerical examples. We 
computed the optimal path numerically for different pa- 
rameters. To find the optimal path one needs to adjust 
a single shooting parameter: the angle between two un- 
stable eigendirections of the endemic fixed point A. Two 
typical examples of numerically computed optimal paths 
are shown in Figs. [1] and [2] [where AK5{1 — SY > 1j ^-iid 
the endemic point is a focus] and in Figs. [3] and |4] [where 
AK5{\ — (5)^ < 1, and the endemic point is a node]. Fig- 
ures [T^ and [3^ show projections of the optimal paths 
on the {x,y) plane. For comparison, they also show the 
mean-field trajectories {px = Py ^ 0) originating in the 
vicinity of the no-disease point x = y = 0, describing an 
epidemic outbreak and approaching the endemic point. 
In contrast to equilibrium systems, the optimal path of 
a large fluctuation is different from the corresponding 
relaxation path. Notice that, although for K = 20 the 
extinction proceeds along a spiral, the difference between 
the two spirals is striking. Figures [T)d and [3)3 show pro- 
jections of the optimal paths on the {px,Py) plane. The 
optimal paths are presented in more detail in Figs. O and 
m where the time dependences of x, y, Px and Py are 
shown. The rescaled action along the optimal path in this 
examples is cr ~ 6.12 x 10"^ for if = 20 and a ~ 0.145 for 
K = 1.8, providing sharp estimates to the logarithm of 
the corresponding mean extinction times of the disease. 

Optimal path and action: asymptotic theory in the 
vicinity of the bifurcation point. For KS <C 1 we can com- 
pute the optimal path and the rescaled action a{K, 5) 
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FIG. 1: (color online) (a) Projection of the optimal path on 
the {x,y) plane (thick black line) and the mean-field trajec- 
tory {px = Py — 0) describing an epidemic outbreak (thin 
red line), (b) Projection of the optimal path on the (px,Py) 
plane, x = S/N - 1, y = I/N; if = 20 and 5 = 0.5. 
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FIG. 2: Optimal path for iC = 20 and 5 = 0.5. Shown are 
X = S/N-1 (a), 1/ = I/N (b), Px=ps-l (c) and py=pi-l 
(d) vs. rescaled time. 



analytically, by exploiting time scale separation. Let 
us introduce rescaled variables: qi = x/S, q2 = yK/d, 
Pi = Px/S'^, and p2 = Py/S. This rescaling is motivated 
by the values of the coordinates of the fixed points, and 
reflects the important feature that, at 5 <C 1, the fluc- 
tuations in the number of susceptible are much weaker 
than the fluctuations in the number of infected. Neglect- 
ing higher order terms in S, we arrive at the following 
approximate equations of motion: 



Pi 



-qi - q2 , 
Pi - q2P2 , 



92 
P2 



KSq2il + qi 
KS {pi -p2- 



l-2p2 



qiP2) 



(10) 



The fixed points become A ~ [—1,1,0,0], B ~ 
[0,0,0,-1], and C = [0,0,0,0]. For KS < 1 the sub- 
system (9i,Pi) is fast, whereas ((72, P2) is slow. On the 
fast time scale (that is, the time scale in the original, 
dimensional variables) the fast subsystem approaches the 
state qi ~ —q2 and pi ~ q2P2 which then slowly evolves 
according to the equations 

92 ^ if5<?2(l-'?2+2p2) , P2 ^ KSp2i2q2-l-p2) (H) 

that are Hamiltonian, as they follow from the reduced 
Hamiltonian Hr{q2,P2) = K5 92^2(1 - 92 +^2)- This 
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FIG. 3: (color online) Same as in Fig. [T] but for K = 1.8. 




FIG. 4; Same as in Fig. [2]but for K = 1.8. 



Hamiltonian appears in the theory of a class of sinqle- 
species models in the vicinity of a bifurcation point . 
As Hr{q2,P2) is independent of time, it is an integral of 
motion. The optimal extinction path goes along the zero- 
energy trajectory 1 — 92+^2 =0 [l^. Evaluating the 
action ([5]) along this line, we find in the leading order: 
So [Nd^/{Kd)] J°P2dq2 = NS^/{2K). For the mean 
extinction time of the disease we obtain 

\n{r)/N ^ Sy{2K) = [m/(2/3)] (1 - r/(3f ; (12) 

this asymptote is valid when iSq ^ 1. 

Dykman et al. [l3| have recently shown that reduced 
Hamiltonian dynamics of the same type as Eqs. (fTT|) 
holds, close to the bifurcation point, in the endemic SIS 
model: still another two-population stochastic epidemic 
model where the infected individuals again become sus- 
ceptible upon recovery. 

In summary, we have developed the eikonal theory for 
stochastic multi-population birth-death systems. The 
theory is especially suitable for analysis of large fluctu- 
ations, such as disease extinction. For the SI and SIR 
models we have found the optimal path to extinction of 
the disease and the mean extinction time. The optimal 
path to extinction, including its remarkable oscillatory 
behavior, is not model-specific. It should be observable 
in stochastic simulations of a broad class of models, and 
in real data on fade out, of infectious diseases in small 
communities. 
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